/**
 * 计算数组最大值
 * @returns {number} 数组中的最大值
 */
Array.prototype.max = function () {
    return Math.max.apply(null, this);
};
/**
 * 计算数组最小值
 * @returns {number} 数组中的最小值
 */
Array.prototype.min = function () {
    return Math.min.apply(null, this);
};
/**
 * 计算数组平均值
 * @returns {number} 数组元素的平均值
 */
Array.prototype.mean = function () {
    let i, sum;
    for (i = 0, sum = 0; i < this.length; i++)
        sum += this[i];
    return sum / this.length;
};
/**
 * 判断点是否在多边形内
 * @param {number} x - 点的x坐标
 * @param {number} y - 点的y坐标
 * @returns {boolean} 如果点在多边形内返回true，否则返回false
 */
Array.prototype.pip = function (x, y) {
    let i, j, c = false;
    for (i = 0, j = this.length - 1; i < this.length; j = i++) {
        if ((this[i][1] > y) !== (this[j][1] > y) &&
            x < ((this[j][0] - this[i][0]) * (y - this[i][1])) / (this[j][1] - this[i][1]) + this[i][0])
            c = !c;
    }
    return c;
};
/**
 * 变差函数模型类
 */
export class VariogramClass {
    /**
     * 构造函数
     * @param {Array<number>} t - 目标值数组
     * @param {Array<number>} x - x坐标数组
     * @param {Array<number>} y - y坐标数组
     * @param {string} model - 变差函数模型类型
     */
    constructor(t, x, y, model) {
        Object.defineProperty(this, "t", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "x", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "y", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "nugget", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "range", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "sill", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "A", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "n", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "K", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "M", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        Object.defineProperty(this, "model", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: void 0
        });
        this.t = t;
        this.x = x;
        this.y = y;
        this.nugget = 0.0;
        this.range = 0.0;
        this.sill = 0.0;
        this.A = 1 / 3;
        this.n = 0;
        this.K = [];
        this.M = [];
        this.model = model;
    }
    /**
     * 高斯变差函数模型
     * @param {number} h - 距离
     * @param {number} nugget - 块金值
     * @param {number} range - 变程
     * @param {number} sill - 基台值
     * @param {number} A - 模型参数
     * @returns {number} 计算得到的变差值
     */
    krigingVariogramGaussian(h, nugget, range, sill, A) {
        return nugget + ((sill - nugget) / range) * (1.0 - Math.exp(-(1.0 / A) * Math.pow(h / range, 2)));
    }
    /**
     * 指数变差函数模型
     * @param {number} h - 距离
     * @param {number} nugget - 块金值
     * @param {number} range - 变程
     * @param {number} sill - 基台值
     * @param {number} A - 模型参数
     * @returns {number} 计算得到的变差值
     */
    krigingVariogramExponential(h, nugget, range, sill, A) {
        return nugget + ((sill - nugget) / range) * (1.0 - Math.exp(-(1.0 / A) * (h / range)));
    }
    /**
     * 球状变差函数模型
     * @param {number} h - 距离
     * @param {number} nugget - 块金值
     * @param {number} range - 变程
     * @param {number} sill - 基台值
     * @returns {number} 计算得到的变差值
     */
    krigingVariogramSpherical(h, nugget, range, sill) {
        if (h > range)
            return nugget + (sill - nugget) / range;
        return nugget + ((sill - nugget) / range) * (1.5 * (h / range) - 0.5 * Math.pow(h / range, 3));
    }
    /**
     * 根据模型类型选择变差函数
     * @param {number} h - 距离
     * @param {number} nugget - 块金值
     * @param {number} range - 变程
     * @param {number} sill - 基台值
     * @param {number} A - 模型参数
     * @returns {number} 计算得到的变差值
     */
    judgeType(h, nugget, range, sill, A) {
        switch (this.model) {
            case 'gaussian': return this.krigingVariogramGaussian(h, nugget, range, sill, A);
            case 'exponential': return this.krigingVariogramExponential(h, nugget, range, sill, A);
            case 'spherical': return this.krigingVariogramExponential(h, nugget, range, sill, A);
            default: return 0;
        }
    }
}
/**
 * 克里金插值主类
 */
export class Kriging {
    constructor() {
        /**
         * 创建填充指定值的数组
         * @param {number} value - 填充值
         * @param {number} n - 数组长度
         * @returns {Array<number>} 创建的数组
         */
        Object.defineProperty(this, "createArrayWithValues", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (value, n) => {
                const array = [];
                for (let i = 0; i < n; i++)
                    array.push(value);
                return array;
            }
        });
        /**
         * 创建对角矩阵
         * @param {number} c - 对角线上的值
         * @param {number} n - 矩阵维度
         * @returns {Array<number>} 对角矩阵
         */
        Object.defineProperty(this, "krigingMatrixDiag", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (c, n) => {
                const Z = this.createArrayWithValues(0, n * n);
                for (let i = 0; i < n; i++)
                    Z[i * n + i] = c;
                return Z;
            }
        });
        /**
         * 矩阵转置
         * @param {Array<number>} X - 输入矩阵
         * @param {number} n - 行数
         * @param {number} m - 列数
         * @returns {Array<number>} 转置后的矩阵
         */
        Object.defineProperty(this, "krigingMatrixTranspose", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (X, n, m) => {
                let i, j;
                const Z = Array(m * n);
                for (i = 0; i < n; i++)
                    for (j = 0; j < m; j++)
                        Z[j * n + i] = X[i * m + j];
                return Z;
            }
        });
        /**
         * 矩阵缩放
         * @param {Array<number>} X - 输入矩阵
         * @param {number} c - 缩放因子
         * @param {number} n - 行数
         * @param {number} m - 列数
         */
        Object.defineProperty(this, "krigingMatrixScale", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (X, c, n, m) => {
                let i, j;
                for (i = 0; i < n; i++)
                    for (j = 0; j < m; j++)
                        X[i * m + j] *= c;
            }
        });
        /**
         * 矩阵加法
         * @param {Array<number>} X - 矩阵1
         * @param {Array<number>} Y - 矩阵2
         * @param {number} n - 行数
         * @param {number} m - 列数
         * @returns {Array<number>} 相加后的矩阵
         */
        Object.defineProperty(this, "krigingMatrixAdd", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (X, Y, n, m) => {
                let i, j;
                const Z = Array(n * m);
                for (i = 0; i < n; i++) {
                    for (j = 0; j < m; j++)
                        Z[i * m + j] = X[i * m + j] + Y[i * m + j];
                }
                return Z;
            }
        });
        /**
         * 朴素矩阵乘法
         * @param {Array<number>} X - 矩阵1
         * @param {Array<number>} Y - 矩阵2
         * @param {number} n - 矩阵1行数
         * @param {number} m - 矩阵1列数/矩阵2行数
         * @param {number} p - 矩阵2列数
         * @returns {Array<number>} 相乘后的矩阵
         */
        Object.defineProperty(this, "krigingMatrixMultiply", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (X, Y, n, m, p) => {
                let i, j, k;
                const Z = Array(n * p);
                for (i = 0; i < n; i++) {
                    for (j = 0; j < p; j++) {
                        Z[i * p + j] = 0;
                        for (k = 0; k < m; k++)
                            Z[i * p + j] += X[i * m + k] * Y[k * p + j];
                    }
                }
                return Z;
            }
        });
        /**
         * 乔里斯基分解
         * @param {Array<number>} X - 输入矩阵
         * @param {number} n - 矩阵维度
         * @returns {boolean} 分解是否成功
         */
        Object.defineProperty(this, "krigingMatrixChol", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (X, n) => {
                let i, j, k;
                const p = Array(n);
                for (i = 0; i < n; i++)
                    p[i] = X[i * n + i];
                for (i = 0; i < n; i++) {
                    for (j = 0; j < i; j++)
                        p[i] -= X[i * n + j] * X[i * n + j];
                    if (p[i] <= 0)
                        return false;
                    p[i] = Math.sqrt(p[i]);
                    for (j = i + 1; j < n; j++) {
                        for (k = 0; k < i; k++)
                            X[j * n + i] -= X[j * n + k] * X[i * n + k];
                        X[j * n + i] /= p[i];
                    }
                }
                for (i = 0; i < n; i++)
                    X[i * n + i] = p[i];
                return true;
            }
        });
        /**
         * 反演乔里斯基分解
         * @param {Array<number>} X - 输入矩阵
         * @param {number} n - 矩阵维度
         */
        Object.defineProperty(this, "krigingMatrixChol2invl", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (X, n) => {
                let i, j, k, sum;
                for (i = 0; i < n; i++) {
                    X[i * n + i] = 1 / X[i * n + i];
                    for (j = i + 1; j < n; j++) {
                        sum = 0;
                        for (k = i; k < j; k++)
                            sum -= X[j * n + k] * X[k * n + i];
                        X[j * n + i] = sum / X[j * n + j];
                    }
                }
                for (i = 0; i < n; i++)
                    for (j = i + 1; j < n; j++)
                        X[i * n + j] = 0;
                for (i = 0; i < n; i++) {
                    X[i * n + i] *= X[i * n + i];
                    for (k = i + 1; k < n; k++)
                        X[i * n + i] += X[k * n + i] * X[k * n + i];
                    for (j = i + 1; j < n; j++) {
                        for (k = j; k < n; k++)
                            X[i * n + j] += X[k * n + i] * X[k * n + j];
                    }
                }
                for (i = 0; i < n; i++)
                    for (j = 0; j < i; j++)
                        X[i * n + j] = X[j * n + i];
            }
        });
        /**
         * 反演高斯-若尔当消元法
         * @param {Array<number>} X - 输入矩阵
         * @param {number} n - 矩阵维度
         * @returns {boolean} 反演是否成功
         */
        Object.defineProperty(this, "krigingMatrixSolve", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (X, n) => {
                const m = n;
                const b = Array(n * n);
                const indxc = Array(n);
                const indxr = Array(n);
                const ipiv = Array(n);
                let i = 0, icol = 0, irow = 0, j = 0, k = 0, l = 0, ll = 0;
                let big, dum, pivinv, temp;
                for (i = 0; i < n; i++) {
                    for (j = 0; j < n; j++) {
                        if (i === j)
                            b[i * n + j] = 1;
                        else
                            b[i * n + j] = 0;
                    }
                }
                for (j = 0; j < n; j++)
                    ipiv[j] = 0;
                for (i = 0; i < n; i++) {
                    big = 0;
                    for (j = 0; j < n; j++) {
                        if (ipiv[j] !== 1) {
                            for (k = 0; k < n; k++) {
                                if (ipiv[k] === 0) {
                                    if (Math.abs(X[j * n + k]) >= big) {
                                        big = Math.abs(X[j * n + k]);
                                        irow = j;
                                        icol = k;
                                    }
                                }
                            }
                        }
                    }
                    ++ipiv[icol];
                    if (irow !== icol) {
                        for (l = 0; l < n; l++) {
                            temp = X[irow * n + l];
                            X[irow * n + l] = X[icol * n + l];
                            X[icol * n + l] = temp;
                        }
                        for (l = 0; l < m; l++) {
                            temp = b[irow * n + l];
                            b[irow * n + l] = b[icol * n + l];
                            b[icol * n + l] = temp;
                        }
                    }
                    indxr[i] = irow;
                    indxc[i] = icol;
                    if (X[icol * n + icol] === 0)
                        return false; // Singular
                    pivinv = 1 / X[icol * n + icol];
                    X[icol * n + icol] = 1;
                    for (l = 0; l < n; l++)
                        X[icol * n + l] *= pivinv;
                    for (l = 0; l < m; l++)
                        b[icol * n + l] *= pivinv;
                    for (ll = 0; ll < n; ll++) {
                        if (ll !== icol) {
                            dum = X[ll * n + icol];
                            X[ll * n + icol] = 0;
                            for (l = 0; l < n; l++)
                                X[ll * n + l] -= X[icol * n + l] * dum;
                            for (l = 0; l < m; l++)
                                b[ll * n + l] -= b[icol * n + l] * dum;
                        }
                    }
                }
                for (l = n - 1; l >= 0; l--) {
                    if (indxr[l] !== indxc[l]) {
                        for (k = 0; k < n; k++) {
                            temp = X[k * n + indxr[l]];
                            X[k * n + indxr[l]] = X[k * n + indxc[l]];
                            X[k * n + indxc[l]] = temp;
                        }
                    }
                }
                return true;
            }
        });
        /**
         * 模型预测
         * @param {number} x - 预测点x坐标
         * @param {number} y - 预测点y坐标
         * @param {VariogramClass} variogram - 变差函数模型
         * @returns {number} 预测值
         */
        Object.defineProperty(this, "predict", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (x, y, variogram) => {
                let i;
                const k = Array(variogram.n);
                for (i = 0; i < variogram.n; i++) {
                    k[i] = variogram.judgeType(Math.pow(Math.pow(x - variogram.x[i], 2) + Math.pow(y - variogram.y[i], 2), 0.5), variogram.nugget, variogram.range, variogram.sill, variogram.A);
                }
                return this.krigingMatrixMultiply(k, variogram.M, 1, variogram.n, 1)[0];
            }
        });
        /**
         * 计算预测方差
         * @param {number} x - 预测点x坐标
         * @param {number} y - 预测点y坐标
         * @param {VariogramClass} variogram - 变差函数模型
         * @returns {number} 预测方差
         */
        Object.defineProperty(this, "variance", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (x, y, variogram) => {
                let i;
                const k = Array(variogram.n);
                for (i = 0; i < variogram.n; i++) {
                    k[i] = variogram.judgeType(Math.pow(Math.pow(x - variogram.x[i], 2) + Math.pow(y - variogram.y[i], 2), 0.5), variogram.nugget, variogram.range, variogram.sill, variogram.A);
                }
                return (variogram.judgeType(0, variogram.nugget, variogram.range, variogram.sill, variogram.A) +
                    this.krigingMatrixMultiply(this.krigingMatrixMultiply(k, variogram.K, 1, variogram.n, variogram.n), k, 1, variogram.n, 1)[0]);
            }
        });
        /**
         * 生成网格矩阵
         * @param {Array<Array<Array<number>>>} polygons - 多边形数组
         * @param {VariogramClass} variogram - 变差函数模型
         * @param {number} width - 网格宽度
         * @returns {IGrid} 网格数据
         */
        Object.defineProperty(this, "grid", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (polygons, variogram, width) => {
                let i, j, k;
                const n = polygons.length;
                if (n === 0)
                    return;
                // 多边形空间边界
                const xlim = [polygons[0][0][0], polygons[0][0][0]];
                const ylim = [polygons[0][0][1], polygons[0][0][1]];
                for (i = 0; i < n; i++) {
                    for (j = 0; j < polygons[i].length; j++) {
                        if (polygons[i][j][0] < xlim[0])
                            xlim[0] = polygons[i][j][0];
                        if (polygons[i][j][0] > xlim[1])
                            xlim[1] = polygons[i][j][0];
                        if (polygons[i][j][1] < ylim[0])
                            ylim[0] = polygons[i][j][1];
                        if (polygons[i][j][1] > ylim[1])
                            ylim[1] = polygons[i][j][1];
                    }
                }
                // 分配O(n^2)空间
                let xtarget, ytarget;
                const a = Array(2);
                const b = Array(2);
                const lxlim = Array(2);
                const lylim = Array(2);
                const x = Math.ceil((xlim[1] - xlim[0]) / width);
                const y = Math.ceil((ylim[1] - ylim[0]) / width);
                const A = {
                    list: Array(x + 1),
                    xlim: Array(2),
                    ylim: Array(2),
                    zlim: Array(2),
                    width: 0
                };
                for (i = 0; i <= x; i++)
                    A.list[i] = Array(y + 1);
                for (i = 0; i < n; i++) {
                    // 多边形[i]的范围
                    lxlim[0] = polygons[i][0][0];
                    lxlim[1] = lxlim[0];
                    lylim[0] = polygons[i][0][1];
                    lylim[1] = lylim[0];
                    for (j = 1; j < polygons[i].length; j++) {
                        if (polygons[i][j][0] < lxlim[0])
                            lxlim[0] = polygons[i][j][0];
                        if (polygons[i][j][0] > lxlim[1])
                            lxlim[1] = polygons[i][j][0];
                        if (polygons[i][j][1] < lylim[0])
                            lylim[0] = polygons[i][j][1];
                        if (polygons[i][j][1] > lylim[1])
                            lylim[1] = polygons[i][j][1];
                    }
                    // 遍历多边形子空间
                    a[0] = Math.floor((lxlim[0] - ((lxlim[0] - xlim[0]) % width) - xlim[0]) / width);
                    a[1] = Math.ceil((lxlim[1] - ((lxlim[1] - xlim[1]) % width) - xlim[0]) / width);
                    b[0] = Math.floor((lylim[0] - ((lylim[0] - ylim[0]) % width) - ylim[0]) / width);
                    b[1] = Math.ceil((lylim[1] - ((lylim[1] - ylim[1]) % width) - ylim[0]) / width);
                    for (j = a[0]; j <= a[1]; j++) {
                        for (k = b[0]; k <= b[1]; k++) {
                            xtarget = xlim[0] + j * width;
                            ytarget = ylim[0] + k * width;
                            if (polygons[i].pip(xtarget, ytarget)) {
                                A.list[j][k] = this.predict(xtarget, ytarget, variogram);
                            }
                        }
                    }
                }
                A.xlim = xlim;
                A.ylim = ylim;
                A.zlim = [variogram.t.min(), variogram.t.max()];
                A.width = width;
                return A;
            }
        });
        /**
         * 在DOM上绘制结果
         * @param {HTMLCanvasElement} canvas - 画布元素
         * @param {IGrid} grid - 网格数据
         * @param {Array<number>} xlim - x轴范围
         * @param {Array<number>} ylim - y轴范围
         * @param {Array<IColor>} colors - 颜色配置
         */
        Object.defineProperty(this, "plot", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (canvas, grid, xlim, ylim, colors) => {
                const ctx = canvas.getContext('2d');
                ctx && ctx.clearRect(0, 0, canvas.width, canvas.height);
                const range = [
                    xlim[1] - xlim[0],
                    ylim[1] - ylim[0],
                    grid.zlim[1] - grid.zlim[0]
                ];
                let i, j, x, y, z;
                const n = grid.list.length;
                const m = grid.list[0].length;
                const wx = Math.ceil((grid.width * canvas.width) / (xlim[1] - xlim[0]));
                const wy = Math.ceil((grid.width * canvas.height) / (ylim[1] - ylim[0]));
                for (i = 0; i < n; i++) {
                    for (j = 0; j < m; j++) {
                        if (grid.list[i][j] === undefined)
                            continue;
                        x =
                            (canvas.width * (i * grid.width + grid.xlim[0] - xlim[0])) / range[0];
                        y =
                            canvas.height *
                                (1 - (j * grid.width + grid.ylim[0] - ylim[0]) / range[1]);
                        z = (grid.list[i][j] - grid.zlim[0]) / range[2];
                        if (z < 0.0)
                            z = 0.0;
                        if (z > 1.0)
                            z = 1.0;
                        if (ctx) {
                            ctx.fillStyle = this.getColor(colors, grid.list[i][j]);
                        }
                        ctx &&
                            ctx.fillRect(Math.round(x - wx / 2), Math.round(y - wy / 2), wx, wy);
                    }
                }
            }
        });
        /**
         * 根据值获取对应颜色
         * @param {Array<IColor>} colors - 颜色配置
         * @param {number} z - 值
         * @returns {string} 颜色字符串
         */
        Object.defineProperty(this, "getColor", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (colors, z) => {
                const l = colors.length;
                for (let i = 0; i < l; i++) {
                    if (z >= colors[i].min && z < colors[i].max)
                        return colors[i].color;
                }
                if (z < 0) {
                    return colors[0].color;
                }
                else {
                    return '';
                }
            }
        });
        /**
        * 使用高斯过程和贝叶斯先验进行训练
        * @param {Array<number>} t - 目标值数组
        * @param {Array<number>} x - x坐标数组
        * @param {Array<number>} y - y坐标数组
        * @param {string} model - 变差函数模型类型
        * @param {number} sigma2 - 噪声方差
        * @param {number} alpha - 正则化参数
        * @returns {VariogramClass} 训练得到的变差函数模型
        */
        Object.defineProperty(this, "train", {
            enumerable: true,
            configurable: true,
            writable: true,
            value: (t, x, y, model, sigma2, alpha) => {
                const variogram = new VariogramClass(t, x, y, model);
                // 计算滞后距离和半方差
                let i, j, k, l, n = t.length;
                const distance = Array((n * n - n) / 2);
                for (i = 0, k = 0; i < n; i++) {
                    for (j = 0; j < i; j++, k++) {
                        distance[k] = [
                            Math.pow(Math.pow(x[i] - x[j], 2) + Math.pow(y[i] - y[j], 2), 0.5),
                            Math.abs(t[i] - t[j])
                        ];
                    }
                }
                distance.sort((a, b) => a[0] - b[0]);
                variogram.range = distance[(n * n - n) / 2 - 1][0];
                // 分箱滞后距离
                const lags = (n * n - n) / 2 > 30 ? 30 : (n * n - n) / 2;
                const tolerance = variogram.range / lags;
                const lag = this.createArrayWithValues(0, lags);
                const semi = this.createArrayWithValues(0, lags);
                if (lags < 30) {
                    for (l = 0; l < lags; l++) {
                        lag[l] = distance[l][0];
                        semi[l] = distance[l][1];
                    }
                }
                else {
                    for (i = 0, j = 0, k = 0, l = 0; i < lags && j < (n * n - n) / 2; i++, k = 0) {
                        while (distance[j][0] <= (i + 1) * tolerance) {
                            lag[l] += distance[j][0];
                            semi[l] += distance[j][1];
                            j++;
                            k++;
                            if (j >= (n * n - n) / 2)
                                break;
                        }
                        if (k > 0) {
                            lag[l] /= k;
                            semi[l] /= k;
                            l++;
                        }
                    }
                    if (l < 2)
                        return variogram; // 错误：点数不足
                }
                // 特征变换
                n = l;
                variogram.range = lag[n - 1] - lag[0];
                const X = this.createArrayWithValues(1, 2 * n);
                const Y = Array(n);
                const A = variogram.A;
                for (i = 0; i < n; i++) {
                    switch (model) {
                        case 'gaussian':
                            X[i * 2 + 1] = 1.0 - Math.exp(-(1.0 / A) * Math.pow(lag[i] / variogram.range, 2));
                            break;
                        case 'exponential':
                            X[i * 2 + 1] = 1.0 - Math.exp((-(1.0 / A) * lag[i]) / variogram.range);
                            break;
                        case 'spherical':
                            X[i * 2 + 1] = 1.5 * (lag[i] / variogram.range) - 0.5 * Math.pow(lag[i] / variogram.range, 3);
                            break;
                    }
                    Y[i] = semi[i];
                }
                // 最小二乘法
                const Xt = this.krigingMatrixTranspose(X, n, 2);
                let Z = this.krigingMatrixMultiply(Xt, X, 2, n, 2);
                Z = this.krigingMatrixAdd(Z, this.krigingMatrixDiag(1 / alpha, 2), 2, 2);
                const cloneZ = Z.slice(0);
                if (this.krigingMatrixChol(Z, 2))
                    this.krigingMatrixChol2invl(Z, 2);
                else {
                    this.krigingMatrixSolve(cloneZ, 2);
                    Z = cloneZ;
                }
                const W = this.krigingMatrixMultiply(this.krigingMatrixMultiply(Z, Xt, 2, 2, n), Y, 2, n, 1);
                // 变差函数参数
                variogram.nugget = W[0];
                variogram.sill = W[1] * variogram.range + variogram.nugget;
                variogram.n = x.length;
                // 带先验的Gram矩阵
                n = x.length;
                let K = Array(n * n);
                for (i = 0; i < n; i++) {
                    for (j = 0; j < i; j++) {
                        K[i * n + j] = variogram.judgeType(Math.pow(Math.pow(x[i] - x[j], 2) + Math.pow(y[i] - y[j], 2), 0.5), variogram.nugget, variogram.range, variogram.sill, variogram.A);
                        K[j * n + i] = K[i * n + j];
                    }
                    K[i * n + i] = variogram.judgeType(0, variogram.nugget, variogram.range, variogram.sill, variogram.A);
                }
                // 投影到目标向量的惩罚Gram矩阵的逆
                let C = this.krigingMatrixAdd(K, this.krigingMatrixDiag(sigma2, n), n, n);
                const cloneC = C.slice(0);
                if (this.krigingMatrixChol(C, n))
                    this.krigingMatrixChol2invl(C, n);
                else {
                    this.krigingMatrixSolve(cloneC, n);
                    C = cloneC;
                }
                // 复制未投影的逆矩阵作为K
                K = C.slice(0);
                const M = this.krigingMatrixMultiply(C, t, n, n, 1);
                variogram.K = K;
                variogram.M = M;
                return variogram;
            }
        });
    }
}
